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ABSTRACT 

We present results from a study of bias and its evolution for galaxy-size halos in a large, high-resolution 
simulation of a low-density, cold dark matter model with a cosmological constant. We consider the 
evolution of bias estimated using three different statistics: two-point correlation function b^, power 
spectrum bp, and a direct correlation of smoothed halo and matter overdensity fields bg. We present 
accurate estimates of the evolution of the matter power spectrum probed deep into the stable clustering 
regime (k ~ [0.1 — 200]/i Mpc -1 at z = 0) and find that its shape and evolution can be well described, with 
only a minor modification, by the fitting formula of Peacock & Dodds (1996). The halo power spectrum 
evolves much slower than the power spectrum of matter and has a different shape which indicates that 
the bias is time- and scale-dependent. At z = 0, the halo power spectrum is anti-biased (bp < 1) with 
respect to the matter power spectrum at wavenumbers k ~ [0.15 — 30]/i Mpc -1 , and provides an excellent 
match to the power spectrum of the APM galaxies at all probed k. In particular, both the halo and 
matter power spectra show an inflection at k rs 0.15/i Mpc -1 , which corresponds to the present-day scale 
of non-linearity and nicely matches the inflection observed in the APM power spectrum. We complement 
the power spectrum analysis with a direct estimate of bias using smoothed halo and matter overdensity 
fields and show that the evolution observed in the simulation in linear and mildly non-linear regimes can 
be well described by the analytical model of Mo & White (1996), if the distinction between formation 
redshift of halos and observation epoch is introduced into the model. We present arguments and evidence 
that at higher overdensities, the evolution of bias is significantly affected by dynamical friction and tidal 
stripping operating on the satellite halos in high-density regions of clusters and groups; we attribute the 
strong anti-bias observed in the halo correlation function and power spectrum to these effects. 

The results of this study show that despite the apparent complexity, the origin and evolution of bias 
can be understood in terms of the processes that drive the formation and evolution of dark matter halos. 
These processes conspire to produce a halo distribution quite different from the overall distribution of 
matter, yet remarkably similar to the observed distribution of galaxies. 

Subject headings: cosmology: theory - large-scale structure of universe - methods: numerical 



1. INTRODUCTION 

The distribution of galaxies may, in general, be biased 
with respect to the overall matter distribution. Therefore, 
the galaxy density field can be used as a probe of matter 
distribution only if we fully understand how the galaxy 
distribution relates to the distribution of matter. Under- 
standing this relationship, the bias, and its evolution are 
of primary importance for the interpretation of the ever 
increasing amount of high-quality galaxy clustering data 
at low and high redshifts. Although the problem of bias 
has been studied extensively during the last decade (e.g., 
Kaiser 1984; Davis et al. 1985; Bardeen et al. 1986; Dekel 
& Silk 1986; Cole & Kaiser 1989; Babul & White 1991), 
new data on galaxy clustering at high redshifts (e.g., at 
z <; 1, Le Fevre et al. 1996; Shepherd et al. 1997; Carl- 
berg et al. 1997; Connolly et al. 1998; Postman et al. 
1998; Carlberg et al. 1998; and, at z ~ 3, Steidel et al. 
1998; Giavalisco et al. 1998; Adelberger et al. 1998) and 
the anticipation of upcoming accurate measurements of 
galaxy clustering at z « from large redshift surveys (e.g., 
Tegmark et al. 1998, and references therein) have recently 
generated substantial theoretical progress in modelling of 
galaxy clustering and bias. 

In hierarchical structure formation models, galaxies are 



assumed to form inside dark matter (DM) halos via the 
energy dissipation by baryons (see Somerville & Primack 
1998 for a recent review). Mo & White (1996) showed how 
bias of dark matter can be predicted analytically in the 
framework of the extended Press-Schechter theory (Bond 
et al. 1991; Bower 1991; Lacey k Cole 1993; see §2). This 
analytical model is rapidly gaining popularity in theoreti- 
cal analyses (see, e.g., recent studies by Kauffman, Nusser 
& Steinmetz 1997, Moscardini et al. 1998, and Baugh et al. 
1998) which requires it to be tested and its capabilities and 
limitations evaluated. The model has been tested against 
numerical simulations by Mo, Jing & White (1996), Cate- 
lan, Matarrese & Porciani (1998), Jing (1998), and Por- 
ciani, Catelan & Lacey (1998), and we present additional 
tests in this paper. More elaborate analytical models have 
been developed by Catelan et al. (1998ab), Porciani et al. 
(1998a), and Sheth & Lemson (1998). 

The effects of non-linearity of the bias on the observable 
statistics have been studied by Fry & Gaztahaga (1993), 
Coles (1993), and Mann, Peacock & Heavens (1997) using 
heuristic models of local non-linear bias. Recently, Coles 
(1993) and Dekel & Lahav (1998) have developed a formal- 
ism for studies of galaxy biasing that allows one to account 
explicitly for the non-linearity and stochasticity of the 
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bias. They have also analyzed the effects of non-linearity 
and stochasticity on results of some of the observational 
analyses. Sherrer & Weinberg (1998) have analyzed effects 
of stochasticity of the local bias on the correlation func- 
tion and power spectrum and concluded that stochasticity 
should not affect the shape of these statistics in the linear 
regime (or, in other words, that in linear regime the bias 
should be linear). Most recently, Narayanan, Berlind & 
Weinberg (1998) studied effects of non-linearity, stochas- 
ticity, and non-locality of the bias using heuristic models 
applied to large iV-body simulations. From the observa- 
tional perspective, Pen (1998) showed how the stochastic- 
ity of the galaxy bias can be tested and measured using 
redshift-space distorsions, and Tegmark & Bromley (1998) 
presented evidence that present-day galaxy bias is non- 
linear and stochastic based on analysis of galaxy clustering 
in the Las Campanas Redshift Survey. 

The evolution of bias in the linear regime has been re- 
cently analyzed by Tegmark & Peebles (1998), who gen- 
eralized the results of Fry (1996) for the case of stochastic 
bias in an arbitrary cosmology. Taruya, Koyama & Soda 
(1998) and Taruya & Soda (1998) used perturbative anal- 
ysis to extend the linear analysis of bias evolution into 
the weakly non-linear regime. Evolution of halo clustering 
and bias in the non-linear regime, has been analyzed in 
several recent numerical studies employing dissipationless 
simulations (e.g., Braincrd & Villumscn 1994; Bagla 1998; 
Colin et al. 1998; Ma 1999), simulations that include both 
dissipationless dark matter and dissipative baryonic com- 
ponents (e.g., Katz, Hernquist & Weinberg 1998; Blanton 
et al. 1998; Cen & Ostriker 1998; Jenkins et al. 1998b), 
and "hybrid" studies in which dissipationless simulations 
are complemented with a semi-analytical model of galaxy 
formation (Roukema et al. 1997; Kauffman et al. 1998ab; 
Diafcrio et al. 1998; Benson et al. 1998; Baugh et al. 
1998; Kolatt et al. 1998). All of these studies, though 
very diverse in their methods, qualitatively agree on one 
important result: the galaxy bias is expected to be non- 
linear, to depend on the properties of the DM halos and the 
galaxies they host, and to be a (generally non-monotonic) 
function of cosmic time. 

In this paper we present results of a detailed analysis of 
halo clustering evolution in a large high-resolution dissipa- 
tionless simulation of a representative and fairly successful 
variant of the cold dark matter (CDM) models: the low- 
density spatially flat CDM model with the cosmological 
constant (ACDM) . The primary goal of this analysis is to 
identify and study the main processes that drive the evolu- 
tion of halo bias in linear and non-linear regimes. Under- 
standing what causes the bias and particular features of 
its evolution (notably, the anti-bias at late stages of evolu- 
tion in some of the models) is crucial for the interpretation 
of current observational and theoretical results. We focus 
therefore on the interpretation of general features of bias 
evolution observed in our simulations and in studies done 
by other authors. The main novel feature of this study 
is inclusion of satellite halos located inside virial radii of 
more massive isolated halos in the halo catalogs; we refer 
reader to Colin et al. (1998) for a more detailed descrip- 
tion of our approach. 

The approach that we have adopted in this project is 
to consider bias estimated using three different statistics: 
two-point correlation function £(r), power spectrum P(k) 



(Peebles 1980), and direct comparison of the smoothed 
halo and matter overdensity (5) fields: 



(1) 



b P (k) = y/Ph{k)/P m {k), (2) 

b s = 5g/5*, (3) 

where quantities with subscripts h and m correspond to 
statistics of the halo and matter distributions, respectively, 
and superscript R indicates the density fields smoothed 
on a scale R. The three estimates of bias given above 
are, of course, related. In the special case of determinis- 
tic local linear bias: = bp = bs- Nevertheless, in the 
general case of stochastic non-linear bias, these estimates 
are complementary to each other and we have chosen to 
consider all three of them to illustrate the manifestations 
and properties of the bias. All three functions, b%, bp, and 
bg, may depend on a number of (both local and non-local) 
parameters; the most important point to notice, however, 
is that they are different functions of their parameters, 
and we use the subscripts to indicate this explicitly. For 
example, b^ and bp have different scale-dependence be- 
cause £(r) and P(k) are related via the Fourier transform, 
£(r) cx J P(k)e l]lL ' T d 3 k, which gives: 



/ b P (k)P m (k)e lk - r d 3 k 
J P m (k)e ik *d 3 k ' 



(4) 



The bias 6f(r) that is scale-dependent in a narrow range 
of r will be scale-dependent in a wide range of k (see Coles 
1993 for a more detailed discussion). We will show below 
that the anti-bias required at small r for the ACDM model 
to be consistent with the z = galaxy correlation func- 
tion, is also perfectly consistent with bp(k) < 1 at small 
wavenumbers (down to k ~ 0.2ft, Mpc -1 ) required for the 
model to be consistent with the galaxy power spectrum 
(Gaztanaga & Baugh 1998; Hoyle et al. 1998). 

The paper is organized as follows. In §2, we give a brief 
account of analytical model of the bias developed by Mo 
& White (1996) and how the epochs of halo formation and 
observation can be separated in this model (e.g., Moscar- 
dini et al. 1998; Catelan et al. 1998a). In §3 and §4 we de- 
scribe the numerical simulation and halo identification and 
selection algorithms used in our analysis. We complement 
the analysis of the evolution of the halo two-point corre- 
lation function based on this simulation and presented in 
Colin et al. (1998) with the analysis of ethe volution of 
the halo power spectrum and bias bp in §5.1. In §5.2 and 
§5.3 we present the analysis of the evolution of bs, esti- 
mated using smoothed halo and matter density fields at 
different epochs, and identify the processes that drive this 
evolution. We discuss the results in §6 and summarize our 
conclusions and arguments in §7. 

2. ANALYTICAL MODEL OF BIAS 

An overdensity field <5(x) = [p(x) — p]/p can be quanti- 
fied in terms of its power spectrum P(k) — (|<5k| 2 ), where 
(5k is the Fourier transform of <5(x). Similarly, if the over- 
density field is smoothed on scale R with a spherically 
symmetric filter W(r, R), the smoothed field can be char- 
acterized by its variance 



a 2 (R) = ({5(x,R)f) = 



(2tt) 3 



P{k)W{R) I d i k, (5) 
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where W is the Fourier transform of the window function. 
In the following and throughout the paper we use the real 
space top-hat mtcr: W(r,R) = (4irR s /Sy^-QiR-r), where 
Q(R—r) is the step function. For this filter the scale R can 
be interchanged with the mass contained within radius R: 
M = (4ir J '3) pR? , where p is the mean density of the uni- 
verse. The standard Press-Schechter model (PS; Press & 
Schcchter 1974) assumes that any region of initial comov- 
ing size R (or mass M) becomes a part of a virialized halo 
by redshift Zf, if its overdensity extrapolated linearly to 
the present epoch is greater than 5 c /D+(zf), where D + (z) 
is the linear growth factor (e.g., Peebles 1980) normalized 
to unity at the present epoch. The value of 5 C is motivated 
by the top- hat collapse model; we use S c = 1.69 throughout 
this work. If the initial overdensity field is gaussian, the 
PS model leads to the following expression for the number 
density of collapsed halos of mass M at redshift z: 



n(M,z,z f )dM 



1 p S c (z,Zf) 



cxp 



2^Ma 3 (M,z) 
$c(z,z f ) 2 
2a 2 (M,z) 



da 2 (M, z) 



dM 
dM, 



(6) 



where 5 c (z,Zf) = S c D + (z)/ D + (zf and a(M,z) is the 
variance of the initial density field smoothed on scale M 
and extrapolated linearly to the epoch z. For the rea- 
sons that will be discussed below, we follow Catelan et al. 
(1998a) in distinguishing the epoch Zf from the "observa- 
tion" epoch z (z < Zf). The Zf- and z-dependencies are 
shown explicitly in the above expression for the mass func- 
tion. Note, however, that n(M, z, Zf) does not change with 
z: for a given power spectrum, the mass function depends 
only on halo mass M and Zf. Equation (6) translates into 
the commonly used form if we assume Zf = z. 

In the seminal paper, Mo & White (1996; hereafter 
MW) showed how the extended Press-Schechter formal- 
ism (EPS; Bond et al. 1991; Bower 1991; Lacey & Cole 
1993) can be used to derive the analytical expression for 
the bias of DM halos. The EPS can be used to derive 
the expression for the conditional mass function of halos 
(Bond et al. 1991). Namely, the number density of ha- 
los of mass M that collapse at epoch zj < z in a region 
of initial Lagrangian radius i?o (mass Mo) and the initial 
overdensity in the growing mode extrapolated linearly to 
the present Sq (So(z) = 5qD + (z)) is given by 



n(M, z,Zf\M ,S )dM 



S 



exp 



2^M[a 2 
[S. 



a 2 } 3 / 2 



da 2 



dM 



where 6 C = 5 c (z,Zf), S = S (z), a 2 = a 2 (M,z), and 
<7q = (T 2 (Mo, z), as defined above. The average overdensity 
of halos of mass M in spheres of overdensity <5o and radius 
Rq at epoch z, 5h{M, z\Ro, Sq), can be obtained by dividing 
number densities given by eqs. (7) and (6) and subtracting 
unity. The halo bias is then defined as b = Sh/So. So far, 
the bias is defined in terms of the Lagrangian radius and 
linearly extrapolated overdensity. For practical purposes, 
however, we need the expression for bias in spheres of ob- 
served radius R and (generally non-linear) overdensity 6. 
MW use a spherical collapse model to relate (R ,S ) to 



(R, S): Rq = (1 + 5)R 3 and Sq are calculated for given R, 
Ro, and d using the equations of spherical collapse. Hav- 
ing made the translation from (R,5) to (i?o,<$o), we can 
calculate the average overdensity of DM halos in spheres 
of the observed radius R and overdensity 5 as 



S h (M,z,z f \R,S) = (1 + 5) 



n(M,z,Zf\Ro,S ) 
n(M, z, zj) 



1, 



b NL (M, z, z f , S) = S h (M, z, z f \R, S)/S, 



(8) 



(9) 



where b^L is the halo bias. Note that in general 6jvl 
depends on the overdensity of matter 5 and is therefore 
non-linear. The quantities (R, S) and (i?o,^o) m eq- (8) 
are related by the spherical collapse model. In the limit of 
linear overdensities and large scales, 5q(z) <C 8 c (z, zf) and 
M ~> M, the bias is given by 



b L (M,z,z f ) = 1 + 



1 



5 c {z,z f y 



(10) 



where v = 5 c (z,Zf)/a(M,z). Eq. (10) shows that in this 
regime the bias is linear (does not depend on 5). Note 
that the bias of halos with a range of masses [M min , M max ] 
should be computed as a mass function weighted average 



b(z,Zf)=n 1 (z,zf) J b(M,z,Zf) n(M,z,Zf) dM, 



(11) 

where n(z, Zf) = f M ma * n(M, z, zf) dM and n(M, z, zf) is 
given by eq. (6). 

We use eqs. (8)-(ll) to calculate the bias predicted by 
this analytical model. The linear overdensity Sq is cal- 
culated for given 8, R, and R = (1 + S)^ 3 R using the 
equations of the spherical collapse model appropriate for 
our ACDM cosmology (e.g., Lahav et al. 1991; Eke et al. 
1996). The resulting function So (5) is well described by 
the fitting formula given by MW (eq.[18] in their paper), 
except for 5 ~ 20 — 30, where deviations reach w 5 — 10%. 

3. NUMERICAL SIMULATION 

We have chosen to study the evolution of the halo bias 
in a representative variant of the CDM-type models: the 
low matter density, flat, CDM model with cosmological 
constant (ACDM): n = 1 - A = 0.3, h = 0.7, where 
O and Oa are present-day matter and vacuum densi- 
ties, and h is the dimensionless Hubble constant defined 
as Ho = lOO/i km/s/Mpc. This model is arguably the 
most successful model in matching a variety of existing 
data. Observations of galaxy cluster evolution (Eke et al. 

1996) , baryon fraction in clusters (Evrard 1997), and high- 
rcdshift supcrnovae (e.g., Pcrlmutter et al. 1998) strongly 
favor the value of matter density fio ~ 0.3, while various 
observational measurements of the Hubble constant (e.g., 
Falco et al. 1997; Salaris & Cassisi 1998; Madore et al. 
1998) tend to converge on the values of h w 0.6 — 0.7. We 
use a normalization of the spectrum of fluctuations that is 
consistent with both observed cluster abundances (Eke et 
al. 1996) and the 4-year COBE data (e.g., Bunn & White 

1997) : ag = 1, where as is the rms fluctuation in spheres 
of Sh^ 1 Mpc comoving radius. 
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The numerical simulation of the ACDM model followed 
the evolution of 256 3 w 1.67 x 10 7 particles in a peri- 
odic 60k' 1 Mpc box. The particle mass is thus w 1.1 x 
10 9 ft _1 M©. The simulation was run using Adaptive Re- 
finement Tree A-body code (ART; Kravtsov, Klypin & 
Khokhlov 1997). The ART code reaches high force res- 
olution by refining all high-density regions with an auto- 
mated refinement algorithm. The refinements are recur- 
sive: the refined regions can also be refined, each subse- 
quent refinement having half of the previous level's cell 
size. This creates an hierarchy of refinement meshes of 
different resolution covering regions of interest. The crite- 
rion for refinement is local overdensity of particles: in the 
simulation presented in this paper the code refined an in- 
dividual cell only if the density of particles (smoothed with 
the cloud-in-cell scheme; Hockney & Eastwood 1981) was 
higher than n t h = 5 particles. Therefore, all regions with 
overdensity higher than S = n t h 2 3L /n, where n is the av- 
erage number density of particles in the cube, were refined 
to the refinement level L. For the simulation presented 
here, n is 1/8. The peak formal dynamic range reached by 
the code on the highest refinement level is 32, 768, which 
corresponds to the smallest grid cell of 1.83ft- -1 kpc; the 
actual force resolution is approximately a factor of two 
lower (see Kravtsov et al. 1997). The simulation that we 
analyze here has been used in Colin et al. (1998), and we 
refer the reader to this paper for further numerical details. 

4. HALO IDENTIFICATION AND SELECTION 

Identification of DM halos in the very high-density envi- 
ronments (e.g., inside groups and clusters) is a challenging 
problem. The goal of this study is to investigate the halo 
bias at both low and high matter overdensitics, and we 
therefore need to identify both isolated halos and satel- 
lite halos orbiting within the virial radii of larger systems. 
The problems associated with halo identification within 
high-density regions are discussed in Klypin et al. (1999; 
KGKK). In this study we use a halo finding algorithm 
called Bound Density Maxima (BDM). The main idea of 
the BDM algorithm is to find positions of local maxima in 
the density field smoothed at a certain scale and to apply 
physically motivated criteria to test whether the identi- 
fied site corresponds to a gravitationally bound halo. The 
detailed description of the algorithm is given in KGKK 
and Colin et al. (1998). The publicly available version 
of the BDM algorithm used here is described in Klypin & 
Holtzman (1997). The halo catalogs used in the present 
study were constructed using numerical parameters given 
in Colin et al. (1998). 

To construct a halo catalog, we have to define selection 
criteria based on particular halo properties. Halo mass is 
usually used to define halo catalogs (e.g., a catalog can be 
constructed by selecting all halos in a given mass range). 
However, the mass and radius are very poorly defined for 
the satellite halos due to tidal stripping which alters a 
halo's mass and physical extent (see KGKK). Therefore, 
we will use maximum circular velocity V max as a proxy 
for the halo mass. This allows us to avoid complications 
related to the mass and radius determination for satellite 
halos. Moreover, when a halo gets stripped V max changes 
less dramatically than the mass, and is therefore a better 
"label" of the halo. For isolated halos, V max and the halo's 
virial mass are directly related. 



TABLE 1 
Halo catalogs 



V max > 120km/s Vmax > 200km/s 



z 


M a . 

vtr 


Nhalo 




^halo 


0.0 


3.0 x 10 11 


4707 


1.4 x 10 12 


1027 


1.0 


2.0 x 10 11 


7867 


9.0 x 10 11 


1443 


2.0 


1.1 x 10 11 


10437 


5.4 x 10 11 


1675 


3.0 


8.0 x 10 10 


9650 


3.5 x 10 11 


1636 



a Masses are given in h 1 Mq. 



For example, a halo with a density distribution de- 
scribed by the Navarro, Frenk & White (hereafter NFW; 
1996, 1997) profile p(r) oc x -1 (l + x)~ 2 (a; = r/R s ; R s is 
the scale-radius): 

2 GM vir c /(2) 

maX Rvir /(c) 2 ' 1 > 

where M V i r and R V i r are the virial mass and radius, 
f(x) = ln(l + x) - x/(l + x), c = R v „/Rs- While for 
the satellite halos V max may not be related to mass in any 
obvious way, it is still the most physically and observa- 
tionally motivated halo quantity. 

We constructed halo catalogs using thresholds in the 
maximum circular velocity (i.e. selecting all halos with 
V max higher than a specified threshold). The cluster-size 
halos are not explicitly excluded from the halo catalogs. 
We assume therefore that the center of each cluster can 
be associated with a central cluster galaxy. The latter 
(due to the lack of hydrodynamics and other relevant pro- 
cesses) cannot be identified in our simulations in any other 
way. We limit the extent of these "galaxies" to the central 
150ft." 1 kpc of the cluster. 

The redshift dependency of the relationship between 
halo mass and V max is expected to evolve because the con- 
centration factor c (see eq. [12]) is expected to evolve with 
redshift. We use the prescription of NFW to compute 
the evolution of c(M, z) and to convert V max to the virial 
mass, but note that this prescription has been found to 
deviate significantly from the results of numerical simula- 
tions (Bullock et al. 1998; Eke, Navarro & Frenk 1998). 
The values of the virial mass corresponding to the V max 
thresholds of 120 and 200 km/s used in our analysis and 
calculated using eq. (12) and the NFW prescription for 
c(M,z) are given in Table 1. This Table also gives the 
number of halos at different epochs identified by the halo 
finder in the simulation box using these thresholds. 

5. RESULTS 

5.1. Evolution of the power spectrum 

Recent numerical and semi-analytical studies have fo- 
cussed on the bias evolution as determined from the 2- 
point correlation function. However, as new, accurate 
measurements of the power spectrum at both low (e.g., 
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Fig. 1. — Evolution of the (a) matter and (b) halo power spectra. Panel (a): the dark matter power spectra, Pm(k), at different 
redshifts (solid lines) are compared with the linear spectra (dotted lines) and predictions of the Peacock & Dodds (1996) fitting 
formula (dashed lines; see text for details). Note that the power spectra in the simulations agree with the analytical predictions 
at all scales, including highly non-linear scales at which the "stable clustering" approximation appears to work well. Panel (b): 
evolution of the power spectrum for halos with maximum circular velocities > 120 km/s (solid lines) as compared to the linear 
evolution of the matter power spectrum (dotted lines). Note that the evolution of the halo power spectrum, Ph(k), is much slower 
than that of the matter spectrum, P m (k), shown in panel (a). At high redshifts the amplitude of Pu(k) exceeds that of P m (k) by 
a factor of ~ 5 — 10, while at lower redshifts the differences are smaller. The ratio of amplitudes, Ph(k) / P m (k) , depends on the 
scale. These differences imply that the halo bias is time- and scale-dependent. 



Baugh & Efstathiou 1993, BE93; Gaztanaga & Baugh 
1998, GB98; Tadros & Efstathiou 1996, TE96; Hoyle et 
al. 1998) and high (Croft et al. 1998; Weinberg et al. 
1998) redshifts become available, it is also useful to ex- 
amine and compare the evolution of the power spectra of 
matter, P m (k), and halo, Ph{k), distributions. 

Figure 1 shows the evolution of real-space P m (k) and 
Ph(k) in our simulation. The halo power spectrum is 
shown for the halo catalog with the maximum circular ve- 
locity threshold of V max > 120 km/s. To estimate the 
power spectra over a wide range of wavenumbers, we have 
used the method of Jenkins et al. (1998a). Both P m (k) 
and Ph(k) have been obtained by combining a series of 
spectra, {Pi(fc)}, in overlapping ranges of k: [k % ~^ x , k % max ], 
where k l max is the maximum wavenumber at which an ac- 
curate estimate of the Pj(fc) can be obtained. To compute 
Pi(k), the computational cube (of size Lf, ox ) is divided 
into i 3 subcubes and the particle (or halo) distributions 
in these subcubes are superposed. The FFT of the result- 
ing density field, estimated using the cloud-in-cell density 
assignment scheme (Hockney & Eastwood 1981), gives an 
accurate estimate of the power spectrum in modes that are 
periodic on scale the L\, ox /i. We have used the FFTs on a 
256 3 grid, and i = 2 m , where m = 0, ...,6 and m = 0, ...,4 
for the dark matter and halos, respectively. Comparisons 
with direct 512 3 -grid FFT spectra suggested the use of 
Kxax = k Ny/ Q i where k l Ny = 128i(2n/L box ) is the Nyquist 
wavenumber for Pj(fe) (k® nax — 2ir/Li, ox ). The individual 



power spectra have been corrected for the shot noise by 
subtracting the noise spectrum estimated using the same 
method. 

Panel (a) of figure 1 shows the evolution of the mat- 
ter power spectrum from z = 5 to the present. For 
comparison we also show the non-linear evolution pre- 
dicted by the fitting formula of Peacock & Dodds (1996, 
hereafter PD96): A 2 NL (k NL ,z) = f NL [A 2 L (k L )], where 
A 2 (fc) = dP(k, z)/d\nk and subscripts L and NL de- 
note linear and non-linear quantities, respectively. The 
analytical expression for /jvl depends on five fitting pa- 
rameters (see §3.3 in PD96) obtained by fitting the power 
spectra of the scale-free iV-body simulations. Although 
there are some minor difficulties in applying these fit- 
ting results to the realistic models with scale-dependent 
power-law spectral index (PD96; Smith et al. 1998; Jenk- 
ins et al. 1998a), the prediction works remarkably well. 
We have been able to match the non-linear spectra in 
our simulation at all epochs with only small change in 
the fitting parameters of PD96. Namely, we have used 
n eff(kL) — d\nP(k)/dhxk\k L (as opposed to an alter- 
native n e ff(ki,/2)) for the estimate of the spectral index 
at wavenumber k^ and slightly changed the power-index 
dependence for the fitting parameter V. This parameter 
controls the amplitude of the high-/c asymptote; instead 
of using V — 11.55(1 + n e f / /3)~ n , with a single time- 
independent value of r) — 0.423 given by PD96, we use rj 
varying from 0.7 at z = 5 to 0.45 at z = 0. The PD96 for- 
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Fig. 2. — Evolution of the halo bias from z = 3.0 to the present. The lower portion of each panel shows linear (dotted lines) and 
non-linear (solid lines) spectra of the matter distribution compared to the halo power spectra for catalogs with two different lower 
cutoffs in maximum circular velocity: > 120 km/s and > 200 km/s (shown by dashed and dot-dashed lines, respectively). The 
upper portion shows the corresponding bias b(k) = [Ph(k) / Pm(k)] 1 ^' 2 . 



mula, with their fitting parameters unchanged, fares well 
at z < 1, but underpredicts the amplitude of the asymp- 
tote at high z. Figure la shows that, with this small 
modification 1 , the PD96 prediction is a success. Never- 
theless, if the desired accuracy of the non-linear spectrum 
estimate is <J 50%, the necessity of such (generally time- 
and cosmology-dependent) modifications should be kept 
in mind. 



Figure lb shows that the evolution of P/ l (fc) is much 
slower than the evolution of the matter power spectrum. 
Although the scale of non-linearity (wavenumber of an up- 
ward inflection of Ph(k)) is seen clearly in the halo spec- 
tra and approximately matches the corresponding scale 
in P m (k) at the same epoch, the shapes of the halo and 
matter power spectra are quite different: the Ph(k), for 
example, can be approximated well by simple power-law, 



This modification is, of course, arbitrary and the same result could possibly be achieved with other changes; for example, by varying the 
wavenumber at which n e f f is computed. 
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while the shape of P m (k) is more complicated. This dif- 
ference, together with the difference in the evolution rate, 
means that the bias of the halo distribution is time- and 
scale-dependent. A similar conclusion has been reached 
by many researchers from comparisons of halo and matter 
two-point correlation functions (e.g., recently, Bagla 1998; 
Colin et al. 1998; Kauffman et al. 1998ab; Katz, Hcrnquist 
& Weinberg 1998; and references therein). Note, however, 
that as we mentioned in §1, the scale-dependence of the 
bias is different in real- and fc-space (see eq.[4]). Although 
we will be referring to both bp(k) and b^(r) (defined in 
eqs. [1] and [2]) as the bias functions, it should be kept in 
mind that they have different functional forms. 

Figure 2 shows the evolution of the halo bias, b(k) = 
[Pft(fc)/P m (fc)] 1 / 2 , from z = 3.0 to the present epoch for 
two halo catalogs selected with low and high-mass thresh- 
olds: V max > 120 km/s and V max > 200 km/s. The bias 
evolves from a value of « 2 — 4 at z = 3 to w 1 at z = 0. 
At high redshifts, the bias depends on the selection thresh- 
old indicating that it is mass-dependent; the differences, 
however, become progressively smaller for lower redshifts. 
Also, the scale-dependence of the bias of the lower-mass 
halos, albeit being redshift-dependent and generally ^ 1, 
is significantly weaker than the scale-dependence of the 
higher-mass halos. At z = 1, for instance, the power 
spectrum of V max > 120 km/s halos follows that of the 
mass almost exactly at all probed wavenumbers. In agree- 
ment with analytical prediction of Scherrer & Weinberg 
(1998), the bias is virtually scale-independent at linear 
scales k < k n i, where k n i is the wavenumber where the 
P m {k) becomes non- linear and exceeds the linear predic- 
tion (k n i w 0.2 — 0.4). Note also that during the evolution 
the bias at these linear scales is driven to the value of 
w 1, as expected in the linear regime (Tegmark & Peebles 
1998). This is true for both low- and high-mass catalogs, 
which indicates that the bias evolution at these scales is 
driven by gravitational growth of clustering that tends to 
erase any initial (mass-dependent) differences in the halo 
and mass distributions (Tegmark & Peebles 1998). 

The bias evolution at non- linear (k > k n i) scales is more 
complicated. The bias evolves to values less than unity 
at z < 1, and its scale-dependence becomes progressively 
stronger over a wider range of wavenumbers. We will dis- 
cuss the possible interpretation of the bias evolution in 
the non-linear regime in the following section. However, 
we would like to point out here that the net result of this 
evolution is the Ph(k) at z = which is significantly anti- 
biased with respect the overall matter distribution but 
which agrees very well with the power spectrum of ob- 
served galaxy distribution. Figure 3 compares the z = 
power spectra of halos and matter in our simulation with 
the power spectrum of galaxies (BE93; GB98; TE96) in the 
APM survey. At small scales (k > 2/iMpc -1 ), the P h (k) 
depends on the catalog's V max threshold. The galaxy 
power spectrum, however, lies comfortably in between 
two likely possibilities of galaxy mass cutoffs. The max- 
imum circular velocities of 120 and 200 km/s correspond 
at z = to the virial masses of « 3 x 10 n ft. -1 M Q and 
ps 1.5 x 10 12 /i -1 Mq, respectively. At larger scales, the 
power spectra of halos and galaxies agree within the errors. 
The "inflection" in the galaxy power spectrum (GB98) is 
reproduced at the correct wavenumber of k s=s 0.2h Mpc -1 . 




k [h Vfpc-'J 

Fig. 3. — Comparison of the z — power spectra of the 
matter (solid line) and halo (dashed and dot-dashed lines) dis- 
tributions with the galaxy power spectra estimated using the 
APM (Gaztanaga & Baugh 1998) and Stromlo-APM surveys 
(Tadros & Efstathiou 1996); the latter is measured in red- 
shift space and must be reduced by ft; 40% to correct out the 
redshift-space effects (see §5.1 for details). The errorbars of 
the galaxy power spectra are 2a. For clarity, the Stromlo- 
APM power spectrum is shown only at k < O.lh Mpc -1 . 
Note that the halo power spectrum, Pf l (k), matches the galaxy 
power spectrum significantly better than does the spectrum of 
the overall matter distribution, P m (k). Both galaxy and halo 
distributions are anti-biased at k ~ 0.2 — 10 Mpc -1 ; the Ph(k) 
at these scales has a shape and amplitude different from those 
of Pm(k). At larger scales, k < 0.2 Mpc -1 , the halo distribu- 
tion is approximately unbiased. 

We interpret the "inflection" wavenumber in the observed 
spectrum as the scale of non-linearity, because both halo 
and matter power spectra in the simulation become non- 
linear at higher k. The amplitude of the matter power 
spectrum at these k, however, is ~ 3 — 4 times higher and 
does not come anywhere close to matching the observed 
power spectrum. This shows clearly that it is crucial 
to consider the distribution of halos, not matter, when 
comparing model predictions to the observations. 

It is not clear whether the power spectrum in the ACDM 
model considered here can reproduce the amplitude and 
shape of the turnover observed in the galaxy power spec- 
trum at k ~ 0.01 - 0.06/i Mpc -1 . Figure 3 shows that 
the ACDM spectrum does not reproduce the turnover in 
the spectrum recovered from the APM angular correla- 
tion function (BE93, GB98); however, it is in better agree- 
ment with the power spectrum derived by TE96 from the 
Stromlo-APM redshift survey. Also, recently published 
power spectrum of the Durham/UKST survey (Hoyle et 
al. 1998) agrees very well with the Stromlo-APM spectrum 
and with the spectrum of the ACDM model. The Stromlo- 
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APM and Durham/UKST spectra have been computed in 
redshift- space, but the differences from the real-space spec- 
trum in the ACDM model at these scales are expected to 
be <; 40% (TE96; Smith ct al. 1998). Even with the 
40% lower amplitude, the spectra are in agreement with 
the model and, within 2a errors, with the APM power 
spectrum. The latter indicates either that there are sys- 
tematic differences between the galaxy surveys or that the 
cosmological model is incorrect (because it predicts an in- 
correct redshift-to-real space correction). There is also 
a possibility that the amplitude and errors of the APM 
power spectrum are somewhat underestimated at these 
large scales (Peacock 1997). Unfortunately, at present, 
statistical and systematic observational errors are too large 
at these scales to be able to make a decisive conclusion. 
In any case, considering the whole range of observationally 
probed wavenumbers, the agreement between the data and 
the model is much better when we compare observations 
to Ph{k), as opposed to P m (k). 

5.2. Overdensity of matter vs. overdensity of halos 

5.2.1. Linear and mildly non-linear overdensities 

The evolution of bias, as determined from the power 
spectra in the previous section, agrees qualitatively with 
the bias evolution derived from the correlation function 
analyses (e.g., recently, Bagla 1998; Colin et al. 1998; 
KaufTman et al. 1998ab; Katz, Hernquist & Weinberg 
1998; and references therein). The bias functions bp(k) = 

^/P h {k)/P m {k) and & c (r) = y/£ h (r)/£ m (r) (see §5.1), de- 
fined using the power spectrum and the correlation func- 
tion, illustrate the scale dependence of the bias. However, 
it is difficult to interpret bp and in terms of the most 
generic definition of bias 2 : Sh = bgS m , where we denote the 
bias in this definition by bg to distinguish it explicitly from 
bp and &£. The bg shows how bias depends on the matter 
density at a fixed scale, the information which cannot be 
extracted easily from b^ and bp. Therefore, to get a full 
picture of the bias evolution, we examine the evolution of 
bg in our simulation. 

To estimate Sh and 5 m , we use the top- hat filter (see §2) 
of comoving radius R = 5/i _1 Mpc. The size of the radius 
is a compromise between halo statistics and the range of 
probed overdensities. Note that our simulation box con- 
tains only 216 independent spheres of this size. This is the 
maximum number of spheres that can be used to study 
the scatter of the bias. However, we are primarily inter- 
ested in the average behavior of bg; therefore, in order to 
probe the wide range of overdensities, we use a large num- 
ber (50, 000) of spheres, placed randomly throughout the 
simulation box. To calculate Sh, we have used halos with 
maximum circular velocities of V ma x > 120 km/s. Due to 
limited mass resolution, the halo catalogs are somewhat 
incomplete for V max <; 100 km/s (see Gottlober, Klypin 
& Kravtsov 1998). The catalog with the threshold of 
120 km/s is complete and contains a large enough num- 
ber of halos (see Table 1) to provide sufficient statistics. 

Figure 4 shows overdensity of dark matter halos, Sh, 
as a function of matter overdensity 5 m at four epochs: 
z = 0, 1, 2, 3. The solid curves show the average relation, 

2 A11 definitions arc, of course, equivalent if the bias is linear. This, however, is not true for the non-linear bias, which appears to be a generic 
feature of the CDM models. 

3 The survival time is the time between halo formation and destruction. 



calculated by averaging Sh of the individual spheres in bins 
of S m . The actual binned distribution of the S m and Sh 
values is shown by grey-scale shades on a grid, where the 
intensity of grey corresponds to the natural logarithm of 
the density of points in the grid cell. The figure shows 
that at all epochs, the halo overdensity is tightly corre- 
lated with the overdensity of matter. However, the slope 
of the correlation, the bias, depends on S m (i.e., the bias 
is non-linear) and changes non-trivially with time. 

The dashed and dot-dashed lines in each panel of Fig- 
ure 4 show predictions of the analytical model of bias de- 
scribed in §2 (namely, eqs. [8] and [9]). To account for 
the range of halo masses used in our halo catalog, we cal- 
culate the mass function weighted bias (eq. [11]) using 
M max = 10 13 /i _1 Mq as an upper limit of integration, and 
rcdshift-dependent M m i„(z) corresponding to our selec- 
tion threshold of maximum circular velocity of 120 km/s 
(see §3 and Table 1). 

The two model curves correspond to different assump- 
tions about formation redshift of halos, Zf. Zf — z (i.e., 
halos forming at the epoch of observation) for the dot- 
dashed curve and Zf = z + 1 for the dashed curve. In 
the standard Prcss-Schechter model, hierarchical build-up 
of halos is a continuous process and, therefore, if mass 
is considered as a defining property of a halo, Zf = z. 
However, this interpretation fails if halos can retain their 
identity over a finite period of time (e.g., see discussion in 
Moscardini et al. 1998). For example, if a halo is accreted 
by a more massive system and orbits in the system's po- 
tential instead of merging instantly, it can be identified at 
z < Zf. The galaxy-conserving model of bias evolution 
(e.g., Moscardini et al. 1998 and references therein) repre- 
sents an extreme in which halos are never destroyed after 
being formed and in which the halo clustering is driven 
solely by the gravitational pull from the surrounding grow- 
ing structures. 

In reality, we expect the merging to take place and the 
halos to have individual merging histories and thus indi- 
vidual survival times 3 (Lacey & Cole 1993; Kitayama & 
Suto 1996). A rigorous treatment of bias should there- 
fore take into account only the halos that survive until the 
epoch of observation z. In practice, this is a difficult task: 
although the halo formation epoch is well defined, defini- 
tion of the halo destruction is not trivial. Lacey & Cole 
(1993) define halo lifetime as the period between formation 
time of a halo and the time by which this halo is incorpo- 
rated into a more massive system. This definition differs 
significantly from how we define the destruction when an- 
alyzing the simulations: the halo is destroyed either when 
it merges with another halo and looses its identity or when 
the tidal stripping brings the mass bound to the halo be- 
low some mass threshold (defined as a selection criterion 
or by the mass resolution of the simulation). Therefore, for 
illustration purposes, we will treat Zf as a free parameter, 
leaving a more rigorous treatment for future work. 

Figure 4 shows that different assumptions about Zf 
lead to significant difference in the behavior of bias pre- 
dicted by eq. (9). Although at S m <; 0.5 the two bias 
predictions are similar, at higher overdensities they di- 
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Fig. 4. — The overdensity of halos Sh vs. the overdensity of matter S m at different epochs (indicated on each panel). Both 
overdensities has been estimated in spheres of radius Rth = 5/i _1 Mpc randomly placed in the simulation box. The solid curves 
show the average relation, calculated by averaging Sh of the individual spheres in bins of 8 m - The actual binned distribution of the 
5 m and Sh values is shown by the shades of grey on a 2D grid, where the intensity of grey corresponds to the natural logarithm of 
the number of spheres in this grid cell. The long-dashed and dot-dashed curves show predictions of the analytical model described 
in §2. The dot-dashed curve is a prediction of the model assuming that formation redshift coincides with observation redshift 
[zf — z), while the long-dashed curve corresponds to assumption Zf = z + 1 (see §5.2 for details). The figure shows that at all 
epochs halo overdensity is tightly correlated with overdensity of matter. However, the slope of the correlation, the bias, depends on 
S m (i.e., the bias is non-linear) and changes non-trivially with time. Note that at z > 1 the halo distribution is biased in overdense 



regions (5 m > 0) but is anti-biased in underdense regions (Sm < 
and almost no bias at low overdensities (see §6 for discussion). 



-0.5). At low-redshift there is an anti-bias at high-overdensities 



verge, the Zf = z + 1 assumption providing a much bet- 
ter match to the simulation results. Note that Zf — z 
underestimates the bias in the simulation at 6 m <; 1. 
This was noted recently by Jing (1998) who studied mass- 
and scale-dependence of bias in the linear regime. For 
the simulation of the ACDM model used in this paper, 
Jing finds that halos of mass 10 n /i -1 M© exhibit bias of 
b 2 (M = lO 11 /*" 1 M Q ) « 0.5 - 0.6, while linear MW bias 
is (eq. [10] assuming z = zi): b\ w 0.28. For the same 
halo mass and Zf = z + 1 = 1, eq. (10) gives b\ w 0.67, 



and therefore the finite survival time of halos may natu- 
rally explain the bias discrepancy. It is worth noting that 
for cluster halos, the br, with Zf = z provides a consider- 
ably better approximation which likely reflects late clus- 
ter formation (z ~ Zf) and/or smaller survival times for 
cluster-size halos, as is indeed predicted by the extended 
Press-Schechter theory (Lacey & Cole 1993). 

5.2.2. Non-linear overdensities 

Figure 5 shows the present-day Sh-S m relation at higher 
overdensities. Although we use a large number of spheres 
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Fig. 5. — The overdensity of halos 8h vs. the overdensity 
of matter S m at z = in spheres of radius -Rth = 5h~ Mpc 
randomly placed in the simulation box. The different sym- 
bols correspond to the spheres in three independent (non- 
overlapping) subcubes of the simulation volume (sub-cube size 
of 30/i." 1 Mpc) that contain most of the massive clusters in the 
simulation. The lines represent the average relations for the 
sub-cubes marked with triangles and crosses. Notice the dif- 
ferences in Sh vs. S m correlation at high S m . We attribute 
these differences to the differences in formation histories of 
structure in these sub-cubes (see §6 for details). The anti- 
bias at very high overdensities (5 m > 10) arises in a region 
surrounding the richest Coma-size cluster in the simulation. 

and oversample the density field, the spheres in indepen- 
dent regions of space are independent and may thus give 
an idea about the true scatter in the halo bias. In order to 
study this scatter, we have divided the computational box 
into eight equal-size (30ft. -1 Mpc) non-overlapping sub- 
cubes. Different symbols (of different colors) in Figure 5 
correspond to spheres in the three sub-cubes that contain 
most of the massive clusters in the simulation at z = 0. 
The other five sub-cubes contain lower density regions and 
are not shown for clarity. 

Figure 5 shows that the differences between 5h-S m re- 
lations in different sub-cubes are significantly larger than 
the scatter of each individual relation. Note, however, that 
even within a single sub-cube there are indications of real 
scatter: the sub-cube represented by the triangle markers 
shows a dichotomy of Sh at fixed 8h ~ 2 — 4. Analysis of 
the sub-cubes has shown that differences in the 5h-5 m cor- 
relation among the sub-cubes are real and are caused by 
the differences in the non-linear structures they contain. 

It appears that the scatter and the differences can be 
explained by the following two effects. The centers of the 
spheres with 6 m ;> 5 fall preferentially in close vicinity to 
the massive group- and cluster-size halos. These cluster 
halos span a wide range of masses and, most importantly, 



formation times. The differences in formation times result 
in different fates of dark matter halos orbiting inside these 
clusters. If, for example, a cluster forms and accretes the 
bulk of its mass and halos early, the halos have time to 
suffer substantial losses of mass due to tidal stripping and 
losses of orbital energy from dynamical friction. Dynam- 
ical friction may lead to a merging between satellite and 
central cluster object, thus resulting in a "loss" of the satel- 
lite. Tidal stripping may lead to a substantial mass loss 
and may ultimately (although at a much slower rate) de- 
crease the halo's maximum circular velocity (see KGKK) 
bringing the halo below our threshold V max . These effects 
are enhanced because the virial radius is smaller at earlier 
epochs and so are the typical pericenters of satellite orbits. 
If, on the other hand, the cluster forms late and accretes 
most of its mass fast and at relatively low redshifts, the 
halos are accreted onto a massive cluster and thus have 
higher orbital energies and orbits with larger pericenters. 
Moreover, many halos simply do not have enough time to 
suffer substantial tidal or orbital energy losses. 

The differences in formation histories may therefore lead 
to significant differences in the halo content among clus- 
ters. This is illustrated in Figure 6 which shows "bias 
profiles" : the ratios of the integral overdensities of halos 
and matter in spheres of increasing radii centered on a 
cluster. The profiles for 5 of the clusters from the three 
sub-cubes of Figure 5 are shown at three different epochs 
z = 0.0, 0.5, 1.0. Note the large differences between profiles 
at z = 1. While the cluster marked CL1 exhibits strong 
anti-bias (i.e., low 6h/5 m ), the halos in cluster CL2 are 
strongly biased. Similar differences are seen in the rest of 
clusters. Interesting differences can also be observed in the 
subsequent evolution of b(r). Cluster CL1 shows very mild 
evolution in b(r) at small (r <J Ih^ 1 Mpc) scales, whereas 
other clusters show very strong evolution between z = 1.0 
and z = 0.5, and much weaker evolution from z — 0.5 until 
present. 

The changes in the rate of evolution may seem coun- 
terintuitive; indeed, the time period corresponding to the 
z = 1.0—0.5 interval is approximately twice as short (~ 2.5 
Gyrs) as the period between z = 0.5 and z = (« 5 Gyrs). 
We could thus expect more significant changes at z < 0.5 
due to more prolonged effects of tidal stripping and dy- 
namical friction. However, as we have noted above, as the 
mass of a cluster grows with time, the halos are accreted 
on higher-pericenter, higher-energy orbits and thus are not 
as likely to approach the dense central region or spend a 
substantial amount of time there. 

To illustrate that this indeed is the case, we have an- 
alyzed the dynamical evolution of halos identified within 
the virial radius of clusters at different moments in time. 
Figure 7 shows the evolution of 5 halos randomly selected 
within the virial radius of cluster CL1 at z ~ 3 (there 
are a total of 10 halos within the virial radius). At all 
epochs, CL1 is the most massive cluster in the simulation 
box; at z = 3 its mass was 1.5 x 10 13 /i _1 M Q . The panels 
in each of the horizontal rows in Figure 7, show the evo- 
lution of particles that bound to the halos at z — 3. The 
orbits of 4 halos are contained within the cluster virial ra- 
dius, s» 150/i _1 kpc, at z = 3 (shown as a circle in each 
panel), where all of these halos merge with the central 
~ lOO/i" 1 kpc size object by z — 1 (i.e., after w 3.5 Gyrs). 
The halo on the most eccentric orbit (bottom row) survives 
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Fig. 6. — Bias profiles (ratio of halo to matter overdensities inside a sphere of radius r centered on the cluster center) of five of 
the rich (M v ir J> lO 14 /)." 1 M©) clusters from the three sub-cubes shown in Fig. 5 at three different epochs z — 0.0,0.5, 1.0. The 
vertical lines indicate the virial radii of clusters at z — 0; M v i T and Z1/2 indicate the z — virial mass and redshift at which cluster 
had half of this mass, respectively. Note the large differences between profiles at z = 1. While cluster marked CL1, the most 
massive cluster in the simulation, exhibits strong anti-bias (i.e., low 8h/S m ), the halos in cluster CL2 are strongly biased. Similar 
differences are seen in the rest of clusters. Interesting differences can also be observed in the subsequent evolution of b(r). Cluster 
CL1 shows very mild evolution in b(r) at small (r < l/i -1 Mpc) scales, whereas other clusters show very strong evolution between 
z = 1 and z — 0.5, and much weaker evolution from z = 0.5 until present. 



until z i=a 0.5 and gets tidally destroyed after that. 

For comparison, Figure 8 shows the evolution of 10 ha- 
los randomly selected within the virial radius of the same 
cluster CL1 at z = 0.5. As before, most of the halos stay 
within the z = 0.5 virial radius (w \ flhr 1 Mpc). However, 
unlike the z — 3 halos, most of them (8 out of 10) survive 
until z = (i.e., during the period of w 5 Gyrs). Although 
some halos suffer substantial mass loss in their outer re- 
gions, the dense halo cores can be identified at z = 0. Note 
that two halos on low-pericenter orbits do merge with the 
central object. 

Finally, we have also visually examined the fate of the 
halos identified in cluster CL5 of Fig. 6 at z = 1 and 
z = 0.5. More than half of the z = 1 halos merge with the 



central object by z — 0.5, while most of the z — 0.5 halos 
survive until z = 0. The difference is caused by both the 
lower typical orbit pericenters at higher rcdshifts and the 
higher efficiency of the dynamical friction due to a smaller 
mass of cluster. It explains the evolution of the bias profile 
shown in Fig. 6. 

To illistrate the relative efficiency of dynamical friction 
at different epochs of cluster evolution, we present the 
evolution of the dynamical friction time-scale in a clus- 
ter. At a given epoch, the dynamical friction time, tf r i c , 
can be estimated using Chandrasekhar's formula (Binncy 
& Tremaine 1987), assuming the cluster mass and density 
distribution and the mass of the satellite. In the right 
column of Figure 9, we present estimates of tf ric for ha- 
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Fig. 7. — The leftmost column of panels shows particles bound to the five halos identified within the virial radius of the cluster 
CL1 shown in Fig. 6 at z — 3 (R v i T ~ 120ft -1 kpc (proper) and M v i r ~ 1.5 x 10 13 /i _1 Mq at z — 3). The rest of the columns show 
positions of the same particles at later moments. Four out of five halos merge with the central cluster halo by z ~ 1. In all panels 
particles are plotted in proper coordinates; the circles in all panels have the same radius equal to the virial radius of the cluster at 

2 = 3. 



los with maximum circular velocity V max = 120 km/s for 
clusters of different final masses. The mass accretion histo- 
ries of clusters of similar present-day mass exhibit a spread 
around an average, typical for this mass, evolution track. 
To account for this spread, we have used both the average 
mass evolution tracks and two individual evolution tracks 
representative of the early- and late-forming tails of the 
population (these tracks represent « 2a deviations from 
the average mass evolution track). The cluster mass evo- 
lution tracks used here (Avila- Reese & Firmani 1997) have 
been generated using the Monte-Carlo method of Lacey & 
Cole (1993). For each epoch, we compute tf ric using eqs. 
(8-10) of KGKK assuming the NFW density distribution 
(with an appropriate c(M,z), see §4) for both the cluster 



and satellite at a distance R V i r /2 from the cluster center ( 
where R V i r is the virial radius of the cluster at this epoch) , 
and accounting explicitly for the mass loss due to the tidal 
stripping. 

The right column of panels in Figure 9 represents the 
dynamical friction time in a different, easier to interpret, 
way. It shows the "merging redshift," defined as the red- 
shift corresponding to the time t + tf r i c , where t is the 
current epoch (redshift z). In the sense of the dynami- 
cal friction time, z merge can be interpreted as a redshift 
at which most of the halos present in cluster at redshift 
2, will merge into a central object. This interpretation 
assumes that the mass of the cluster would not change, 
which is not correct. The actual value should therefore be 
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Fig. 8. — The same as in Fig. 7 but for ten halos randomly selected from the halos identified within the virial radius of the 
cluster at z = 0.5 {R vlr fa 1.1ft -1 Mpc (proper) and M v i r « 5.8 x 10 14 ft -1 Mq). The first and third columns from the left show 
the positions of the particles bound to the halos at z = 0.5, while the second and the fourth panels show the positions of the same 
particles at z = 0. Note that despite substantial mass losses, eight out of ten of these halos can be identified at z — as distinct 
dense clumps of particles. The circles in all panels have the same radius equal to the virial radius of the cluster at z = 0.5. 



considered as an upper limit on the actual z merge . Nega- 
tive values z merger < mean that the dynamical friction 
time is larger than the time span between redshift z and 
the present. 

Figure 9 shows that regardless of the final cluster mass, 
dynamical friction is efficient at z J> 3. For small final mass 
clusters, dynamical friction is efficient throughout the en- 
tire cluster evolution, while for medium- and high-mass 



clusters, dynamical friction effectively switches off as the 
cluster exceeds a certain mass threshold. Cluster CL1 in 
Figure 6, evolves very close to the average mass evolution 
track of the 6 x 10 14 ft _1 Mq final mass clusters. From 
the estimate of z merge we can therefore expect that all 
of its satellite halos will merge into the central object by 
z ~ 0.5 — 1. Note also that for all cluster masses, only halos 
accreted before 2 = 0.5 can be substantially influenced by 
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Fig. 9. — Estimates of the dynamical friction time-scale as a function of time for clusters of different final (z — 0) mass Mo- The 
right column presents estimates of tf r ic for satellite halos of maximum circular velocity V m ax = 120 km/s. The three curves cor- 
respond to the average evolution (solid curves), early formation (dashed curves) and late formation (dot-dashed curves) computed 
using Monte-Carlo realizations of mass accretion histories (see §5.3 for details). The dotted curve shows the age of the Universe in 
the ACDM model studied here as a function of redshift. The left column shows the corresponding evolution of "merging redshift" 
defined as a redshift corresponding to t(z) +tf ric (z), where t(z) is age of the Universe at redshift z and tf ric is the corresponding 
dynamical friction time-scales (both quantities are shown in the right panel). The curve marking has the same meaning as in the 
right panel. Negative values z mer ger < mean that the dynamical friction time is larger than the time span between redshift z and 
the present. 



dynamical friction. However, even at z < 0.5, dynamical 
friction may remain efficient in less massive clusters and 
groups, which, if accreted by a cluster, will be depleted of 
halos and will tend to lower the bias value further in this 
cluster. Both of these results are in qualitative agreement 
with simulation (see Figs. 7,8). 

6. DISCUSSION 

Results and arguments presented in the previous section 
suggest that we can identify major processes that drive 
the evolution of the halo bias. In the linear (S m ^ 1; 
the overdensities quoted here and below are for a density 



field smoothed with the top-hat filter of radius Rth = 
5/i _1 Mpc) and mildly non-linear regimes (S m <; 3) the 
analytical model of bias developed by Mo & White (1996; 
see §2) is in good agreement with the results of our simu- 
lation, if the formation epoch of halos z / is distinguished 
from the epoch of observation z (or, in other words, if halos 
are assumed to retain their identity during a finite inter- 
val of time after their formation). This model reproduces 
the non-linearity and time evolution of the bias observed 
in the simulation well (see Fig. 4). The evolution of the 
bias in linear and quasi-linear regimes is, therefore, driven 
primarily by the halo collapse and merging rates specific 
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for a given cosmological model. 

Interplay between halo formation and the merging rate 
in high-density regions (affecting n(M, z, Zf\Ro, do) in 
eq.[8]) and halo formation and the merging rate in the 
field (affecting n(M,z,Zf)) leads to the decrease of the 
bias with time for halos in a given mass range. This is sim- 
ply because halos of a given mass collapse earlier in high- 
density regions than they do in the field; the high-redshift 
objects therefore represent rare events in the density field 
and are initially strongly clustered due to modulation by 
large-scale modes. The number of halos in high-density 
regions is then decreased due to merging, while number 
density of halos in the field may still increase (or level off 
depending on mass range considered) at lower redshifts. 
The ratio is thus a decreasing function of time. In the re- 
gions of negative overdensity (5 m <^ —0.5), the evolution 
is reverse: formation of halos of a given mass in these re- 
gions is delayed and at early epochs their number density 
is below the average (see Fig. 4). The halos in under- 
dense regions are thus anti-biased at early epochs. Figure 
4 shows that bias at 5 m < —0.5 increases during evolu- 
tion as the collapse threshold is reaching lower and lower 
overdensities and more and more halos are being formed 
in these undcrdcnsc regions (Fig. 5 shows that b w 1 for 
— 1 <J S m <J 2 at z = 0). The prediction of the analytical 
model at these overdensities matches our numerical results 
nicely. 

In high-density regions (S m J> 3), the evolution and 
amplitude of bias appears to be significantly affected by 
dynamical friction 4 . Indeed, the process results in halo 
mergers with the central cluster halo and reduces thereby 
the ratio of halo to matter overdensities. The efficiency 
of dynamical friction depends sensisitively on the mass of 
satellite halos, properties of halo orbits, and mass of the 
cluster. The latter evolves rapidly with time and switches 
dynamical friction off at some epoch (z <~ 0.5 — 1.0 for the 
ACDM model studied here). The mass accretion history 
is a stochastic process and some scatter in the mass evo- 
lution is expected for clusters that have the same mass at 
z = (Lacey & Cole 1993): some clusters accrete most 
of the mass early on, while others accrete most of their 
mass at lower redshifts. Therefore, for a given observation 
epoch z, different clusters may be affected by the dynam- 
ical friction to a different degree. For example, a cluster 
which have accreted most of its mass (and halos) just prior 
to z, will be less affected than a cluster of the same mass 
that accreted most of its halos earlier because dynamical 
friction had more time to operate in this cluster. The sit- 
uation is even more complicated, because clusters, even 
those in which dynamical friction becomes inefficient due 
to mass increase, may accrete smaller clusters and groups 
which, in turn, have different evolution histories and there- 
fore different ratios of halo to matter overdensities. 

Dynamical friction is not the only process that affects 
halo counts in clusters and groups. Tidal fields of clusters 
strip the outer parts of satellite halos which may result in 
substantial mass loss for a medium- and high-mass clus- 
ters (up to a factor of 5 — 20 depending on paramaters of 



the halo orbit and the period of time the halo spends in 
cluster; see, e.g., KGKK). The maximum circular velocity, 
V max , of the halo changes only mildly, even in the case of 
severe mass loss (which, as a reminder, is the reason we 
use it for the halo selection). Nevertheless, for some halos 
V m ax may decrease below the selection threshold, in which 
case these halos will "drop out" of the catalogs. This pro- 
cess likely contributes to the mild evolution of bias seen in 
Fig. 6 at z < 0.5. The efficiency of tidal stripping is lower 
for lower-mass systems; therefore, we expect it to be im- 
portant only in relatively massive (M„j r ;> 10 14 /i _1 M Q ) 
clusters. While this effect may seem to depend on the par- 
ticular halo selection procedure used in our analysis, sim- 
ilar effects may arise for other selection procedures 5 . It is 
clear, for example, that this effect would be even more se- 
vere had we chosen to select halos using their bound mass. 
Observationally, galaxy catalogs are usually selected using 
a fixed luminosity limit in a given wavelength band. If the 
luminosity of galaxies in this band evolves differently in 
clusters than in the field (which is strongly suggested by a 
variety of observations), the uniform selection criterion is 
bound to select somewhat different galaxy populations in 
high- and low-density regions. 

Different evolution histories of different systems may re- 
sult in different locally evaluated bias. If, for example, 
matter and halo density fields are smoothed at some suffi- 
ciently large (larger than typical cluster size) scale, regions 
of a similar matter overdensity may correspond to different 
halo overdensities, because the latter depends on the evo- 
lution history of systems encompassed within the smooth- 
ing scale. We can expect, therefore, that bias evaluated 
at a finite smoothing scale will exhibit some scatter in dif- 
ferent regions of space. Note that this scatter arises not 
from the "stochasiticity" of the halo formation, but from 
the fact that the same matter overdensity may correspond 
to regions of very different evolution histories and, thus, 
of different halo content. On the other hand, the differ- 
ences in the evolution histories result from the modulation 
by large-scale modes in the density distribution. The bias 
is thus also modulated by the large-scale modes and is 
therefore non-local. Note that this "stochasticity" should 
decrease, as one smoothes density fields on progressively 
larger scales because effects of the modulation by large- 
scale modes are smaller on larger scales. In general, if the 
relation between halo and matter density is non-linear, the 
bias estimated from a density field smoothed at any partic- 
ular scale will be non-linear and will exhibit some scatter 
(Dckel & Lahav 1998). In our simulation, differences in 
the bias between different regions of the computational 
volume seen in Fig. 5 are caused by different numbers 
and formation histories of clusters and groups in these re- 
gions. The regions that exhibit weaker bias (i.e., stronger 
anti-bias) are the regions that contain clusters with earlier 
formation epochs. The Coma-size cluster that already had 
mass of w 1.3 x 10 13 /i _1 M Q at z = 3 exhibits the highest 
matter overdensity and the strongest anti-bias. 

Clusters and groups of galaxies contribute a signifi- 
cant fraction of the galaxy clustering signal at small, 



4 It is clear that dynamical friction is important for halo evolution; indeed, binary halo mergers are also due to dynamical friction. However, 
here we discuss dynamical friction that operates on satellites orbiting inside a more massive system (a group or a cluster), which leads to the 
decay of their orbits, and, ultimately, to a merger with the central cluster object. 

5 Note that we are bound to use some selection procedure because in most cases we are interested in studying the clustering of a particular 
(selected) class of objects: halos of a given type, galaxies of a given luminosity or color, etc. 
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r <^ 1 — 2h Mpc, scales and significantly affect cluster- 
ing amplitude at larger scales. The anti-bias arising due to 
dynamical processes in groups and clusters can therefore 
explain the anti-bias seen in comparisons of halo and mat- 
ter two-point correlation functions and power spectra. Our 
results then imply that understanding of the evolution of 
small-scale galaxy clustering requires a deeper understand- 
ing of the evolution of galaxies in groups and clusters. 

The major caveat in interpretation of these dissipation- 
less results is a correspondence between dark matter halos 
and observable galaxies. We note, however, that the def- 
inition of dark matter halo in this study is significantly 
different from a conventional definition, which eliminates 
many problems of halo-galaxy mapping related to the over- 
merging problem (see KGKK and Colin ct al. 1998 for 
discussion). In the framework of the hierarchical struc- 
ture formation modelled here, it seems likely that in ev- 
ery sufficiently massive (M ^ lO 11 ^" 1 M Q ) gravitationally 
bound halo, at some epoch baryons will cool, form stars, 
and produce an object ressembling a galaxy (e.g., White 
& Rees 1978; Kauffman, Nusser & Stcinmetz 1997; Yepes 
et al. 1997). This is indeed a cornerstone assumption 
of the semi-analytical models of galaxy formation (e.g., 
Somcrville & Primack 1998). Therefore, even though we 
cannot predict detailed properties of galaxies hosted by 
dark matter halos (which would require inclusion of a sub- 
stantial amount of additional physics in the simulations), it 
is likely that the overall halo distribution should be repre- 
sentative of the expected distribution of the overall galaxy 
population. Excellent agreement between clustering of 
galaxy-size halos in our simulation and observed galaxy 
clustering, demonstrated in comparisons of the two-point 
correlation functions (Colin et al. 1998) and power spectra 
(see Fig. 3), is an indirect but strong indication in favor 
of this point. 

Anti-bias, similar in amplitude and scale-dependence to 
that observed in our simulation, has been found in other 
recent numerical studies that employ either non-adiabatic 
hydrodynamics (Jenkins ct al. 1998b) or semi-analytic 
recipes (Kauffman et al. 1998ab; Diaferio et al. 1998; 
Benson et al. 1998) to model galaxy formation and evolu- 
tion of the two-point correlation function. These studies 
are very different in their modelling technqiues and the 
agreement indicates that anti-bias is real and is not a nu- 
merical artefact of a particular simulation. Anti-bias was 
traditionally an unfavored proposition because it is easier 
to envision a biased rather than anti-biased galaxy forma- 
tion. However, as we have argued above, the anti-bias may 
arise naturally during the dynamical evolution of the halo 
population, even though halo formation is biased. Indeed, 
all of the numerical studies mentioned above are similar to 
the present study in explicit modelling of the orbital evolu- 
tion of halos in groups and clusters. Diaferio et al. (1998) 
present bias profiles b(r) of the clusters in their simulation 
that are in good qualitative agreement with the profiles 
shown in Fig. 6, indicating that similar dynamical pro- 
cesses are probably causing the anti-bias in their study. 

We can therefore expect that simulations affected by 
the overmerging problem and in which a recipe for split- 
ting or weighting the overmerged halos is used may over- 
estimate the clustering amplitude at small-scales and not 
show (or show a weaker) anti-bias. The amount of anti- 
bias should also depend on the box size because small-size 



30 — 50ft. -1 Mpc) boxes are unlikely to contain clus- 
ters in the high-mass tail of the mass function. Such a 
trend is indeed observed; the comparison of bias in the 
ACDM simulations of 30/i -1 Mpc and 60/i _1 Mpc boxes 
presented in KGKK and Colm et al. (1998) shows that 
anti-bias is stronger in the 60/i _1 Mpc simulation. Finally 
and most importantly, the sensitivity of the small-scale 
bias amplitude to the abundance and evolution histories 
of clusters indicates that we can expect some differences 
between cosmological models. The models that form clus- 
ters at systematically later epochs and/or in lesser abun- 
dance, should exhibit a weaker anti-bias (or even absence 
of anti-bias), than models that form clusters earlier and in 
larger numbers. Thus, for example, a tCDM model ap- 
pears to result in a higher value of z = bias than the 
ACDM model (Kauffman et al. 1998ab, Colm et al. 1998; 
Diaferio et al. 1998). A more systematic large-box study 
is needed, however, to test and quantify such dependence 
on cosmology. 

We would like to note that the anti-bias (in the amount 
predicted by the numerical studies) is actually needed to 
reconcile the otherwise very successful ACDM model with 
observed galaxy clustering at z = (e.g., Klypin, Pri- 
mack & Holtzman 1996; Cole et al. 1997; Jenkins et al. 
1998a). Note that this requirement applies to the over- 
all galaxy population represented in large galaxy surveys 
such as the CfA and APM; sub-populations of galaxies 
may exhibit significantly different clustering properties, as 
is indicated by the differences between clustering of spiral 
and elliptical galaxies (e.g. Guzzo et al. 1997), infrared- 
and optically-selected galaxies (e.g. Hoyle et al. 1998), 
etc. In the simulation presented here, certain sub-samples 
of halos are clustered more strongly than the overall pop- 
ulation. For example, as shown in Gottlober, Klypin & 
Kravtsov (1998), the population of halos that loose mass 
at z < 1 (the halos that are accreted and being tidally 
stripped by clusters) are actually biased at z — with re- 
spect to the matter distribution, as opposed to the strongly 
anti-biased distribution of the entire halo population. The 
simulation used in our analysis was also used by Kolatt 
et al. (1998), who showed that interpretation of the high 
clustering amplitude of the high-redshift galaxies is not 
unique. The amplitude can be reproduced equally well 
when these galaxies are assumed to be located in high- 
mass halos, or are assumed to have smaller masses but 
undergo a starburst due to a collision/merger with another 
galaxy. Recent studies by Blanton et al. (1998), Cen & 
Ostriker (1998), Kauffmann ct al. (1998ab), and Diaferio 
et al. (1998) demonstrate the existence of age, luminos- 
ity, color segregation of clustering in their models. Some 
differences in clustering properties of certain sub-samples 
from the overall population can therefore be expected. 

Another observational requirement for the ACDM 
model is a scale-dependent bias (e.g., Jenkins et al. 1998a): 
the shape of the observed galaxy correlation function is 
rather different from the shape of the non-linear matter 
correlation function; this scale-dependency is also repro- 
duced nicely in the simulations. Recent studies of galaxy 
power spectrum by Gaztahaga & Baugh (1998) and Hoyle 
et al. (1998) stress that the anti-bias is required at rather 
low wavenumbers. However, as we have noted above (see 
§1) and have illustrated in Fig. 3, the amount of small- 
scale anti-bias observed in the galaxy-correlation function 
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(Colin et al. 1998) is sufficient to produce the required 
anti-bias in the power spectrum at low wavenumbers. In- 
deed, the power spectrum of the halo distribution in our 
simulation matches nicely the power spectrum of the APM 
galaxies at all probed wavenumbers. The matter power 
spectrum, on the other hand, is different from the galaxy 
power spectrum at a very high (> 5 — 10cx) significance 
level. The existence of non-linear and scale-dependent 
bias of the galaxy distribution may affect other analy- 
ses that either assume there is no bias or that the bias 
is linear. These include estimates of the matter density 
SI from peculiar velocities and redshift distorsions (Dekel 
& Lahav 1998; Pen 1998) and from the observed mass 
to light ratios in galaxy groups and clusters (Diaferio et 
al. 1998), estimates of the baryon density in the Uni- 
verse (e.g., Goldberg & Strauss 1998; Meiksin, White & 
Peacock 1998), estimates of the cosmological parameters 
based on the joint analysis of galaxy redshift surveys, cos- 
mic microwave background anisotropies, and high-redshift 
supernovae (e.g, Eisenstein, Hu & Tegmark & 1998), etc. 
In this respect, the lesson of the present analysis is that 
any use of the galaxy density field as a cosmological probe 
requires very careful testing and evaluation. 

7. CONCLUSIONS 

We presented results from a study of the origin and evo- 
lution of bias of the dark matter halo distribution in a 
large, high-resolution simulation using a low-matter den- 
sity, flat, CDM model with cosmological constant. The 
following conclusions can be drawn from the results pre- 
sented in this paper. 

1 . The evolution of the power spectrum of the halo dis- 
tribution is significantly slower than the evolution of the 
matter power spectrum at all (both linear and non-linear) 
scales. The halo and matter power spectra also have signif- 
icantly different shapes. The differences in shape and rate 
of evolution imply time- and scale-dependent bias of the 
halo distribution which is in qualitative agreement with 
the results of the correlation function analyses. We stress, 
however that the scale dependence of the bias determined 
from the power spectrum bp(k) = \J Ph(k) / P rn (k) is dif- 
ferent from the scale dependence of (r) = y / £,h(r)/£,m(r) 1 
because the two statistics are related via the Fourier trans- 
form (see §1). Put simply, bp(k) cannot be obtained from 
b^(r) by a naive substitution of variable r oc l/k. At z = 
the halo power spectrum in our simulation matches the 
observed power spectrum of the APM galaxies well. 

2. Despite the differences in shape, the power spectra of 
both matter and halos exhibit a distinct "inflection point" 
at approximately the same wavenumber, corresponding to 
the scale of non-linearity (i.e., the scale at which the power 
spectra begin to deviate significantly from the linear power 
spectrum). The inflection scale is ~ 0.15 — 0.2h Mpc -1 
and coincides with the inflection observed in the power 
spectrum of APM galaxies (Gaztanaga & Baugh 1998); 
therefore, we interpret the inflection in the APM power 
spectrum as the present-day scale of non-linearity fcjvi- In 
the simulation analyzed here, fcjvi evolves with time from 
- lh Mpc" 1 at z = 5 to w 0.15-0.2/j Mpc" 1 at z = 0. We 
should note that the distinct inflection point can be seen 
only in the real-space power spectrum; the non-linear am- 
plitude of the redshift-space power spectrum is strongly 



suppressed and the inflection in the power spectrum of 
both matter and halos is smoothed out (see Figs. 3 & 4 in 
Gottlober et al. 1998). 

3. The analytic fitting formula of Peacock & Dodds 
(1996), with only minor tuning, provides an excellent 
match to both the shape and evolution rate of the mat- 
ter power spectrum in our simulation. The latter probes 
deep into the non-linear regime, down to wavenumbers of 
<~ 200ft- Mpc" 1 (at z = 0); we find that clustering of dark 
matter in the highly non-linear regime in the simulation is 
approximately stationary (stable clustering). 

4. In addition to bp, we examine the evolution of bias 
defined using smoothed halo and matter overdensities (5h 
and S m ) as bs = Sh/S m . In agreement with results from 
the correlation function and power spectrum analyses, we 
find that bs is non-linear (i.e., depends on 5 m ) and time- 
dependent. If we modify the model and assume that ha- 
los can retain their identity for a finite period of time 
after their formation and distinguish between formation 
and observation epochs Zf and z, the analytic model of 
Mo & White (1996) can describe both the non-linearity 
and evolution of bs at linear and quasi-linear overdensities 
(S m ^5 5 — 7, here and below the overdensities are smoothed 
with the top-hat smoothing radius of 5ft" 1 Mpc). The 
original model (z — zj) significantly underestimates the 
bias of galaxy-size halos at 6 m ^ 1. 

5. We argue that at non-linear overdensities the bias 
evolution is significantly affected by dynamical friction and 
tidal stripping of halos in groups and clusters. Both pro- 
cesses tend to reduce the number density of cluster halos 
of a given mass range, thereby reducing the bias and re- 
sulting in anti-bias at z ^ 0.5 at 8 m J> 5 in the ACDM 
model studied here. The effect of dynamical friction de- 
pends sensitively on the cluster formation history, which 
introduces a certain degree of scatter into the bias bs- 

In summary, the evolution of the bias of galaxy-size ha- 
los observed in the simulation in linear and quasi-linear 
regimes can probably be fully described using the extended 
Press-Schechter theory. In other words, the evolution of 
bias in this regime results from an interplay between halo 
formation and merging rates in different regions and in 
the field. In the non-linear regime, the halo bias evolution 
appears to be driven by the dynamical processes inside 
clusters and groups. Thus, despite the apparent complex- 
ity, we believe that the origin and evolution of bias can be 
understood in terms of the processes that drive the forma- 
tion and evolution of dark matter halos and galaxies that 
they host: collapse from the density peaks, merging, tidal 
stripping and morphological transformation in the high- 
density regions. Our results show that these processes may 
conspire to produce a halo distribution quite different from 
the overall distribution of matter, yet remarkably similar to 
the observed distribution of galaxies. This result implies 
that detailed modeling of the small-scale galaxy cluster- 
ing requires a good understanding of galaxy evolution in 
clusters. We would like to emphasize, therefore, the im- 
portance of further efforts in modeling galaxy evolution in 
both clusters and in the field. 
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